
% 判断是否为外部调用
callparm = exist('pca_l1_calloutside', 'var');
if callparm == 1
    if pca_l1_calloutside == 1
        figswitch = 0;
    end
else
    clear;clc;close all;
    figswitch = 1;
end

addpath('D:\Radar_Simulation\SidelobeLevel\');
%%

avrpow = @(s, n) s'*s/n; 
rcm = @(s, n) sqrt(sum((abs(s)).^6) / n) / sqrt((avrpow(s, n))^3);

Os = 4;
% 载波数量
N = 64;

mask = ones(N, 1);

% % 信号子载波位置随机
infnum = 1;
r1 = randperm(N, infnum);
% load("r1fortest.mat", "r1");
for i = 1:infnum
    mask(r1(i), 1) = 0;
end
qam_map = 1;
carriermask = zeros(sum(mask), 1);
carriermask(1:end) = 1.0; %%——完全子载波覆盖；
% carriermask(1:end) = 0.8; %%——无任何子载波覆盖；
% carriermask(1:end) = 0.6; %%——完全子载波覆盖；
% carriermask(1:end) = 0.4; %%——无任何子载波覆盖；
% carriermask(1:end) = 0.2; %%——完全子载波覆盖；
% carriermask(1:end) = 0.0; %%——无任何子载波覆盖；
% carriermask(65:65) = 40;


% % 信号子载波居中
% infnum = 16;
% mask((N / 2 - infnum / 2 + 1) : (N / 2 + infnum / 2),1) = 0;
% % infnum = N - sum(mask);
% qam_map = 16;
% carriermask = zeros(sum(mask), 1);
% % carriermask(1:end) = 20; %%——完全子载波覆盖；
% % carriermask(1:end) = 0; %%——无任何子载波覆盖；
% carriermask(55:68) = 40;
% % carriermask(1:end) = 1.0; %%——完全子载波覆盖；
% % 以B/2对称分布的起始子载波功率，偏离信息子载波为offset，单边起始子载波宽度wd

% % 信号子载波靠左边
% infnum = 16;
% mask(1:infnum,1) = 0;
% carriermask = zeros(sum(mask), 1);
% % 完全子载波覆盖
% offset = 0;
% wd = 112;
% carriermask(1+offset:offset+wd) = 1;
% % carriermask(1:8) = 1;
% % carriermask(17:24) = 1;
% % carriermask(65:72) = 1;
% D = sum(mask);
% alpha = 0;
% for i = 1:length(carriermask)
% % carriermask(i) = 1 / (1 + alpha * (i - 1)^2 / D^2);
% % carriermask(i) = 0;
% end

% 信号子载波靠右边
% mask(113:128,1) = 0;
% infnum = N - sum(mask);
% carriermask = zeros(sum(mask), 1);
% % % 完全子载波覆盖
% offset = 64;
% wd = 16;
% carriermask(sum(mask)-offset-wd+1:sum(mask)-offset) = 1;

% 信号子载波靠边
% mask(1:8,1) = 0;
% mask(121:128,1) = 0;
% carriermask = zeros(sum(mask), 1);
% offset = 0;
% wd = 56;
% carriermask(1+offset:offset+wd) = 1;
% carriermask(113-offset-wd:112-offset) = 1;


pha = pi * 2 * rand(N, 1);
% load('phafortest.mat', 'pha');


wi = zeros(Os*N, infnum);
fi = zeros(infnum, Os*N);
wr = zeros(Os*N, (N - infnum));
fr = zeros((N - infnum), Os*N);
mapi = zeros(infnum, N);
remapi = zeros(N, infnum);
mapr = zeros((N-infnum), N);
remapr = zeros(N, (N-infnum));

xr_index = (find(mask))';
xi_index = (find(1-mask))';
x_index = (1:N);
outx_index = (N+1:Os*N);

% iindex = 0;
% rindex = 0;
% for i = 1:N
%     if(mask(i, 1) == 0)
%         iindex = iindex + 1;
%         wi(:,iindex) = exp(2i * pi / Os / N * (i - 1) * (0:Os*N-1)') / Os / N;
% %         fi(iindex,:) = exp(-2i * pi / Os / N * (i - 1) * (0:Os*N-1)');
% %         mapi(iindex, i) = 1;
% %         remapi(i, iindex) = 1;
%     else
%         rindex = rindex + 1;
%         wr(:,rindex) = exp(2i * pi / Os / N * (i - 1) * (0:Os*N-1)') / Os / N;
% %         fr(rindex,:) = exp(-2i * pi / Os / N * (i - 1) * (0:Os*N-1)');
% %         mapr(rindex, i) = 1;
% %         remapr(i, rindex) = 1;
%     end
% end

wo = zeros(Os*N, (Os*N - N));
iindex = 1;
for i = (Os*N - N+1):Os*N
    wo(:,iindex) = exp(2i * pi / Os / N * (i - 1) * (0:Os*N-1)') / Os / N;
    iindex = iindex+1;
end

% rweights的产生
rweights = exp(1i * pha(xr_index));
if(norm(rweights, 2) ~= 0)
    rweights = rweights / norm(rweights, 2);
end
rweights = carriermask .* rweights;

rweights = 1 / norm(rweights, 2) * rweights;

% rweights(:,:) = 0;
% 使用QPSK的方式产生iweights
%   在rweights归一化的情况下，iweights功率占比是Numi/N时，iweights的幅度
% amp = sqrt(1 / sum(mask));
% iweights = amp * rand_qpsk(N- sum(mask));

% 使用16-QAM作为调制信息
% iweights = qam(16, infnum);
% iweights = sqrt(infnum / (N-infnum)) / norm(iweights, 2) * iweights;

% 使用8-PSK作为调制信息
% iweights = psk(8, infnum);
% iweights = sqrt(infnum / (N-infnum)) / norm(iweights, 2) * iweights;

load("iweightsfortest.mat", 'iweights');
% iweights(:,:) = 0;

iweightsamp = abs(iweights);

% % 参数设置
weights = zeros(N, 1);
weights(xr_index) = rweights;
if ~isempty(xi_index)
    weights(xi_index) = iweights;
end

% load("Init_weightsrandom_N_128_mod_8-PSK_mc_5000_wInit_1.mat", 'exp_init_weights');
% windex = 2;
% weights = (exp_init_weights(2*windex-1,:)).';
% xi_index = exp_init_weights(2*windex, 1:infnum);
% xr_index = exp_init_weights(2*windex, infnum+1:end);
% iweights = weights(xi_index);
% rweights = weights(xr_index);

Os_len = Os * N;
wi = exp(2i *pi / Os_len * (0:Os_len-1)' * (xi_index - 1)) / Os_len;
wr = exp(2i *pi / Os_len * (0:Os_len-1)' * (xr_index - 1)) / Os_len;

allpwr = weights'*weights;
freepwr = weights'*weights - iweights'*iweights;

startweights = weights;
s = ifft(weights, Os*N);
starts = s;
os16StartS = ifft(weights, 16*N);
disp(['cve of init signal is ', num2str(cve(starts))]);

Xr = wr';
Xi = wi';
Ad = Xr' * Xr;
XX = diag(diag(Ad));

if ~isempty(xi_index)
    c = wi * iweights;
else
    c = zeros(Os*N, 1);
end


% b = wr*rweights + c;
% b = c;
% startweights = zeros(N, 1);
% startweights(xr_index) = rweights;
% if ~isempty(xi_index)
%     startweights(xi_index) = iweights;
% end
% angleofstartb = angle(startweights);
% figure(44);
% plot(angleofstartb);
% return

b = starts ./ abs(starts);
startb = b;
scalestartb = startb * sqrt(N / (N - infnum)) / Os / N;
fftofstartb = fft(scalestartb);
powAuxVectorInFreq =  zeros(N,1);
if ~isempty(xi_index)
    powAuxVectorInFreq(xi_index) = fftofstartb(xi_index);
end
powAuxVectorInFreq(xr_index) = fftofstartb(xr_index);
figure(43)
yyaxis left
plot(abs(powAuxVectorInFreq))
yyaxis right;
plot(abs(fftofstartb))
xlim(gca, [0 128])

%% 迭代次数
rounds = 1000;
delta = zeros(rounds,1);
% rweights = X*b / norm(X*b, 2);
% weights = remapi*iweights + remapr*rweights;
% s = ifft(weights, Os*N);


fftdeltabarr = zeros(rounds, (N - infnum));

timeResidualComponentsArr = zeros(rounds, Os*N);
freqResidualComponentArr = zeros(rounds, N - infnum);

% 测试代码-每次地带后的分布--start--
viewstart = 1;
viewend = Os*N;
figure(49);hold off;
plot3((1:Os*N), real(starts), imag(starts));hold on;
% plot3((1:Os*N), real(c), imag(c));
plot3([1;Os*N], [0; 0], [0; 0], 'k', 'lineWidth', 1.5);
% plotx = zeros(1,3*Os*N);
% ploty = zeros(1,3*Os*N);
% plotz = zeros(1,3*Os*N);
% plotx(1:3:end) = (1:Os*N);
% plotx(2:3:end) = (1:Os*N);
% plotx(3:3:end) = (1:Os*N);
% ploty(1:3:end) = 0;
% ploty(2:3:end) = real(s.');
% ploty(3:3:end) = 0;
% plotz(1:3:end) = 0;
% plotz(2:3:end) = imag(s.');
% plotz(3:3:end) = 0;
% 
% plot3(plotx, ploty, plotz, '-b', 'lineWidth', 1.5);
%
view([0.5,-1,-1]);
xlim(gca, [viewstart, viewend]);
ylim(gca, [-1e-3, 1e-3]);
zlim(gca, [-1e-3, 1e-3]);
%     plot(s, '*');
%     axis(1e-4*[-8 8 -8 8]);
% 测试代码-每次地带后的分布--end--

% accelerate method pre-work --start--
rweights = Xr * b;
temp = angle(b);
normcoe = 1 / norm(rweights, 2);
rweights = normcoe * rweights;
orgrweights = rweights;
weights(xr_index) = rweights;
% b = wr * rweights;
% s = b + c;
s = ifft(weights, Os*N);
b = s ./ abs(s);

b = exp(1i * 2 * pi * rand(Os*N, 1));
non_smoth_weights = (fft(b));

orgb = b;
orgAngleOfb = angle(b);
orgphas = angle(b);
newAngleOfb1 = orgAngleOfb - temp;
lastdelta = zeros(N - infnum, 1);
% accelerate method pre-work --end--

%% 门限设置
weightPowTh = 1.44 * (rweights' * rweights) / length(rweights);
weightTh = sqrt(weightPowTh);

% b = c./abs(c);

CM = rcm(s, length(s));
disp(['CM of init s:' num2str(20*log10(CM))]);
Ex4 = sum((abs(s)).^4);
disp(['E(x4) of init s:' num2str(Ex4)]);
disp(['cve of init s:' num2str(10*log10(cve(s)))]);
disp(['power of init s:' num2str(s'*s)]);
disp('----------------------------------------');

if ~isempty(xi_index)
    iweights = 1.0 * iweights;
    weights(xi_index) = iweights;
    c = wi * iweights;
end

angChanges = zeros(rounds, Os*N);
% bmodulus = zeros(rounds, Os*N);
% xrangles = zeros(rounds, N-infnum);
% xrmodulus = zeros(rounds, N-infnum);
xrrecord = zeros(rounds, N-infnum);
srecord = zeros(rounds, Os*N);
brecord = zeros(rounds, Os*N);
cvesrocord = zeros(rounds, 1);

timmeanpow = sqrt(s' * s / Os / N);
deltasum = zeros(size(b));

K = 90;
temphis = zeros(Os*N, K);
shis = complex(zeros(Os*N, K));
deltaangewithchis = zeros(Os*N, rounds);

ifftmask = ifft(mask, Os*N);
rpow = sqrt(1/Os/N);
% ifftmask(10:(Os_len-8)) = 0;

tic;
for i = 1:rounds

    
% % 测试代码-每次迭代的相位分布-start--
%     figure(48);hold on
%     phaseofb = angle(b);
%     plot(phaseofb(1), phaseofb(Os+1), '*')
% % 测试代码-每次迭代的相位分布-end--
    
% %  basic method --start--

%     orgb = b;
%     rweights = X * b;
% %     orgrweights = rweights;
%     normcoe = 1 / norm(rweights, 2);
%     rweights = normcoe * rweights;
%     orgrweights = rweights;
% %     b = wr * rweights;
% %     b = b + c;
% %     b = b ./ abs(b);
% % %     
%     weights(xr_index) = rweights;
%     s = ifft(weights, Os*N);
%     b = s ./ abs(s);
% %     cves = cve(s);
% %     disp("iteration counts " + num2str(i) + ...
% %         ", cve after iteration1:" + num2str(10 * log10(cves)) + "dB/"+num2str(cve(s)));

%  basic method --end--

%  accelerate method
%     rweights = X * b;
%     normcoe = 1 / norm(rweights, 2);
%     rweights = normcoe * rweights;
% %     rweights = rweights + 0.9 * (rweights - orgrweights);
% %     normcoe = 1 / norm(rweights, 2);
% %     rweights = normcoe * rweights; 
% 
% %     lastdelta'*(rweights - orgrweights)/(norm(lastdelta, 2)^2)
% %     lastdelta = rweights - orgrweights;
% 
%     weights(xr_index) = rweights;
%     orgrweights = rweights;
%     s = ifft(weights, Os*N);
%     b = s./abs(s);
%     
% %   accelerate part
%     newphas = angle(b);
%     delphas = newphas - orgphas;
%     b = exp(1i * (newphas + 1 * delphas));
%     orgphas = newphas;
%     
%     cves = cve(s);
%     disp("iteration counts " + num2str(i) + ...
%         ", cve after iteration1:" + num2str(10 * log10(cves)) + "dB/"+num2str(cve(s)));
% %   accelerate method end

% %  max subcarrier weight constraint start
    orgb = b;
    orgAngleOfb1 = angle(b);
    orgAngleOfb = orgAngleOfb1;
    orgs = s;
    rweights = Xr * b;
    normcoe = 1 / norm(rweights, 2);
    rweights = normcoe * rweights;
    iweights = Xi * b;
    
    moduleOfRWeights = abs(rweights);
    tempmax = zeros(size(moduleOfRWeights));
    maxposindex = [];
    temppow = 1;
    
    if(weightTh)
    
        while(1)
    %         [maxweight, maxpos] = max(moduleOfRWeights);
    %         
    %         if(maxweight <= weightTh)
    %             break;
    %         else
    %             maxposindex = [maxposindex, maxpos];
    %             moduleOfRWeights(maxpos) = 0;
    %             temppow = temppow - weightPowTh;
    %             
    %             moduleOfRWeights = sqrt(temppow) / norm(moduleOfRWeights) * moduleOfRWeights;
    %         end

            indices = find(moduleOfRWeights > weightTh);
            if(isempty(indices))
                break;
            end
            temppow = temppow - length(indices) * weightPowTh;
            moduleOfRWeights(indices) = 0;
            tempmax(indices) = weightTh;
            moduleOfRWeights = sqrt(temppow) / norm(moduleOfRWeights) * moduleOfRWeights;

        end
    %     moduleOfRWeights(maxposindex) = weightTh;
        moduleOfRWeights = moduleOfRWeights + tempmax;
        rweights = (moduleOfRWeights ./ abs(rweights)) .* rweights;
    end
    
%     xrangles(i,:) = angle((rweights - orgrweights));
%     xrmodulus(i,:) = abs((rweights - orgrweights));
    xrrecord(i,:) = rweights - orgrweights;
    
    weights(xr_index) = rweights;
%     weights(xi_index) = iweightsamp(1) * (iweights ./ abs(iweights));
    orgrweights = rweights;
    s = ifft(weights, Os*N);
    
% %     PAR方法
%     s = sqrt(Os*N) / norm(s, 2) * s;
%     maxDiv = max(abs(s) - 1);
%     th = 0.5*maxDiv;
%     upth = 1+th;
%     for j = 1:Os*N
%         
%         b(j) = s(j);
%         
%         if(abs(s(j)) > upth)
%             b(j) = upth * s(j) / abs(s(j));
%         end
%         
%     end
    
%     PVR方法

%     s = sqrt(Os*N) / norm(s, 2) * s;
%     maxDiv = max(abs(abs(s) - 1));
%     th = 0.5*maxDiv;
%     upth = 1+th;
%     lowth = 1-th;
%     for j = 1:Os*N
%         
%         b(j) = s(j);
%         
%         if(abs(s(j)) > upth)
%             b(j) = upth * s(j) / abs(s(j));
%         end
%         
%         if(abs(s(j)) < lowth)
%             b(j) = lowth * s(j) / abs(s(j));           
%         end
%     end

%     CM方法的b更新
    b = s./abs(s);

    
    angleOfb = angle(b);
    
    brecord(i,:) = b;
    temp = (angleOfb - orgAngleOfb1);
    
    for k = 1:K-1
        temphis(:,k) = temphis(:,k+1);
        shis(:,k) = shis(:,k+1);
    end
%     temphis(:,K) = temp;
    shis(:,K) = s;

    deltasum = deltasum + temp;
    
    %     if mod(i, 20)==0
    %         deltasum(:,:) = 0;
    %     end
    
    
    
%     figure(41);hold off;
%     plot(wrapToPi(temp) / pi, 'r');hold on;
%     plot(wrapToPi(deltasum) / pi, '--b');
% %     ylim([-0.1 0.1])
%     grid on;
    
    os16s = ifft(weights, 16*N);
    cves = cve(os16s);
    cvesrocord(i,1) = cves;
%     disp("iteration counts " + num2str(i) + ...
%         ", cve after iteration1:" + num2str(10 * log10(cves)) + "dB/"+num2str(cves));
    
% % %     accelerate part
    
    newAngleOfb = angleOfb - orgAngleOfb;
    angChanges(i, :) = newAngleOfb;
    orgAngleOfb = angleOfb;
%     angleOfb = angleOfb + 1.0 * deltasum;
    tempang = orgAngleOfb1 + 1 * deltasum + 1 * temp;
    angleOfb = angleOfb + 1.0 * newAngleOfb;
%     b = exp(1i * angleOfb);


    
% % %     accelerate part 2
    
%     dels = shis(:,K) - shis(:,K-1);
%     s = s + dels;
%     
%     if(i <=50)
%         b = exp(1i * tempang);
%     else
%         b = s ./ abs(s);
%     end
    
% %  max subcarrier weight constraint end

% % 纯时域操作-start
% 
%     tempfilter = cconv(ifftmask, b, Os*N);
%     tempfilter = rpow / norm(tempfilter, 2) * ifftmask;
%     
% %     abstempfilter= abs(tempfilter);
% %     figure(55)
% %     hold on;
% %     plot(abstempfilter)
% %     hold on
% %     plot(abs(c))
% %     xlim([45 50])
% 
% %     tempfilter1 = tempfilter;
% % %     tempfilter1(1) = 0;
% %     tempfilter1 = cconv(tempfilter1, b, Os*N);
% %     tempfilter1 = b .* conj(tempfilter1);
% %     abstempfilter1 = abs(tempfilter1);
% %     angtempfilter1 = wrapToPi(angle(tempfilter1));
% %     c_m_b = b .* conj(c);
% %     deltaangewithc = wrapToPi(angle(c_m_b));
% %     deltaabswithc = abs(c_m_b);
% %     deltaangewithchis(:, i) = deltaangewithc;
% %     figure(55)
% %     yyaxis left;
% %     hold off;
% %     plot(deltaangewithc, 'b-');
% %     hold on;
% %     plot(angtempfilter1, 'r-')
% %     ylim([-pi pi])
% %     yyaxis right;
% %     hold off;
% %     plot(deltaabswithc, 'b--');
% %     hold on;
% %     plot(abstempfilter1, 'r--');
%     
%     tempfilter = cconv(tempfilter, b, Os*N);
%     tempfilter = tempfilter + c;
%     orgb = b;
%     b = tempfilter ./ abs(tempfilter);
%     
%     deltaangewithchis(:, i) = wrapToPi(angle(b .* conj(orgb)));
%     
% 
% 
% 
%     rweights = X * b;
% %     normcoe = 1 / norm(rweights, 2);
% %     rweights = normcoe * rweights;
%     weights(xr_index) = rweights;
%     s = ifft(weights, Os*N);


% % 纯时域操作-end
    
    % 测试代码-每次地带后的分布--start--
%     figure(49);hold off;
%     plot3((1:Os*N), real(s), imag(s));hold on;
%     plot3((1:Os*N), real(c), imag(c));hold on;
%     plot3([1, Os*N], [0, 0], [0, 0], 'k', 'lineWidth', 1.5);hold on;
%     temp = sqrt(s'*s / Os / N);
%     temps = temp * s ./ abs(s);
%     if i == 1
%         tempsstart = temps;
%         starttemps = temp * startb ./ abs(startb);
%         plot3((1:Os*N), real(starttemps), imag(starttemps), 'g');hold on;
%     end
%     plot3((1:Os*N), real(temps), imag(temps));hold on;
%     tempdeltas = s - temps;
% %     plot3((1:Os*N), real(tempdeltas), imag(tempdeltas), 'm');
%     view([0.5,-1,-1]);
%     xlim(gca, [viewstart, viewend]);
%     ylim(gca, [-1e-3, 1e-3]);
%     zlim(gca, [-1e-3, 1e-3]);
    
%     plotx = zeros(1,3*Os*N);
%     ploty = zeros(1,3*Os*N);
%     plotz = zeros(1,3*Os*N);
%     plotx(1:3:end) = (1:Os*N);
%     plotx(2:3:end) = (1:Os*N);
%     plotx(3:3:end) = (1:Os*N);
%     ploty(1:3:end) = 0;
%     ploty(2:3:end) = real(s.');
%     ploty(3:3:end) = 0;
%     plotz(1:3:end) = 0;
%     plotz(2:3:end) = imag(s.');
%     plotz(3:3:end) = 0;
% 
%     plot3(plotx, ploty, plotz, '-b', 'lineWidth', 1.5);
% 
%     view([0.5,1,1]);
%     plot(s, '*');
%     axis(1e-4*[-8 8 -8 8]);
    % 测试代码-每次地带后的分布--end--
    
    % 测试代码-Bussagang分解--start--
%     tempout = b;
%     tempin = orgb;
% %     tempin = unwrap(angle(orgb));
% %     tempout = unwrap(angle(b));
%     timeResidualComponentsArr(i, :) = tempout - tempin'*tempout/(norm(tempin, 2)^2)*tempin;
%     tempin = orgrweights;
%     tempout = rweights;
%     freqResidualComponentArr(i, :) = tempout - tempin'*tempout/(norm(tempin, 2)^2)*tempin;
%     freqResidualComponentArr(i, :) = tempin'*tempout/(norm(tempin, 2)^2);
    % 测试代码-Bussagang分解--end--
    
    % 测试代码--start--
%     temp1 = fft(b);
%     temp2 = fft(orgb)-temp1;
%     prj = (temp1(outx_index))'*(temp2(outx_index))/(norm(temp2(outx_index), 2)^2);
%     disp("Projection is : " + num2str(prj));
%     figure(47);hold off;
%     plot(abs(temp1));hold on;
%     plot(abs(prj*temp2));
%     plot(abs(prj*temp2 - temp1));
    % 测试代码--end--
    
    % 使用相关性，获得算法加速--start--
%     fftdeltab = orgrweights - newrweights;
%     fftdeltabarr(i, :) = fftdeltab;
%     rweights = newrweights - 2*fftdeltab;
%     normcoe = 1 / norm(rweights, 2);
%     rweights = normcoe * rweights;
%     weights(xr_index) = rweights;
%     s = ifft(weights, Os*N);
%     cves = cve(s);
%     disp("iteration counts " + num2str(i) + ...
%         ", cve after iteration2:"+num2str(10 * log10(cves))+"dB/"+num2str(cve(s)));
%     b = s./abs(s);
    % 使用相关性，获得算法加速--end--
%     disp(" ");

% wb = conj(b) .* (Xr'*Xr*b + norm(Xr*b)*Xi'*iweights);
% figure(81);
% % plot(abs(wb));
% stem(real(wb), '-r');hold on;
% stem(imag(wb), '--b');hold off;

end

tempweisig = ifft(weights);
tempweisig = tempweisig ./ abs(tempweisig);
smoth_weights = fft(tempweisig);

wb = conj(b) .* (Xr'*Xr*b + norm(Xr*b)*Xi'*iweights);
figure(81);hold on;
% plot(angle(wb));
stem(abs(wb), ':m', 'LineWidth', 2);
stem(imag(wb), 'r');
stem(real(wb), '--b');
wrn = norm(Xr(:,1),2)^2;
stem((1:length(wb)), wrn*ones(1, length(wb)), ':g', 'LineWidth', 2);

% load('starbwithpeakpsd.mat', 'startb');

h = b .* conj(startb);
rhh = (h * h').';

% figure(61)
% hold off
% temp = zeros(1, Os_len);
% for i=1:Os_len
%     temp = temp + (circshift(rhh(i, :), (1-i)));
%     plot(real(temp));
%     if(mod(i, 100) == 0)
%         100==100;
%     end
% end

% figure(62)
% hold off
% angles = wrapToPi(angle(rhh));
% 
% for i = 1:Os_len
%     angles = (circshift(angles, (i - 1), 2));
%     plot(diag(angles));
%     ylim([-pi pi]);
% end

figure(63)
hold off

for i = 1:rounds
    histogram(deltaangewithchis(:,i), 32);
%     xlim([-pi pi]);
    title("第"+num2str(i)+"次迭代各元素旋转角度的分布");
end

ffth = ifft(h);
figure(63)
plot(abs(ffth))
ffthh = ffth * ffth';


temp = temp.' / Os_len;
ph = fft(temp);
fftstartb = fft(startb);
pbin = fftstartb .* conj(fftstartb);
fftoutb = fft(b);
pbout = fftoutb .* conj(fftoutb);
figure(59)
hold on;
plot(pbout, 'b-');
temp2 = cconv(pbin, ph, Os_len) / (Os_len);
plot(abs(temp2), 'r-');
plot(abs(pbin), 'g-.');
plot(abs(0.7945*pbin), 'm-.');
figure(54)
plot(abs(fftstartb));
hold on;
plot(abs(fftoutb));
xlim([0 256])

runtime = toc;
avragetime = runtime / rounds;

% figure(42)
% for i = 1:N*Os
%     plot(unwrap(angChanges(:, i) / pi));
% end
% plot(cvesrocord);

% angRelMet = angChanges * angChanges';
% for i = 1:rounds
%     angNorms = norm(angChanges(i,:), 2);
%     angRelMet(i,:) = angRelMet(i,:) / angNorms;
%     angRelMet(:,i) = angRelMet(:,i) / angNorms;
% end


% 测试代码-Bussagang分解--start--
% residualCorr = timeResidualComponentsArr*timeResidualComponentsArr';
% residualCorr = freqResidualComponentArr*freqResidualComponentArr';
% for i=1:rounds
%     residualCorr(i,:) = residualCorr(i,:) / (norm(freqResidualComponentArr(i,:), 2)^2);
%     residualCorr(i, 1:i-1) = 0;
% end
% 测试代码-Bussagang分解--end--

% return
% rweights = X*b / norm(X*b, 2);
% weights(xr_index) = rweights;
os16OutS = ifft(weights, 16*N);
os16Finalcve = cve(os16OutS);
avr = norm(os16OutS, 2) / sqrt(16*N);
enve = abs(os16OutS);
convunitrat = avr ./ enve;
disp(['oversample enhanced cve is ', num2str(10 * log10(os16Finalcve)), '/', num2str(os16Finalcve)]);

acf_start = conv(os16StartS, flip(conj(os16StartS)));
acf_out = conv(os16OutS, flip(conj(os16OutS)));
[~, psll_start] = psll(os16StartS);
[~, psll_out] = psll(os16OutS);
isll_out = isll(os16OutS);
disp(['ACF of init/out wave, PSLL: ' num2str(psll_start) '/'...
    num2str(psll_out) '(dB), ISLL: ' num2str(isll_out) '.'])


start_powspec = (abs(startweights)).^2;
out_powspec = (abs(weights)).^2;

del_powspec = mapr * (out_powspec - start_powspec);
% del_powspec = abs(del_powspec);
del_powspec = del_powspec ./ (mapr * start_powspec);
% 记录频谱改变程度
del_ampspec = mapr * ((abs(weights) - abs(startweights)));
del_ratio_ampspec = mapr * ((abs(weights) - abs(startweights)) ./ abs(startweights));
del_pha = (angle(b .* conj(startb)));
cum_start_powspec = cumsum(start_powspec);
cum_out_powspec = cumsum(out_powspec);

% % 研究两次迭代之间，xr上变化量的相关关系
delxrrela = zeros(rounds, rounds);
for i = 1:rounds
    for j = i:rounds
        delxrrela(i,j) = conj(fftdeltabarr(i,:)) * (fftdeltabarr(j,:)).' / (norm(fftdeltabarr(i,:), 2)^2);
    end
end

CM = rcm(s, length(s));
disp(['CM of result s:' num2str(20*log10(CM))]);
PVR = max(abs(s)) / min(abs(s));
disp(['PVR of result s:' num2str(PVR)]);
PAR = max(abs(s)) / norm(s, 2) * sqrt(Os*N);
disp(['PAR of result s:' num2str(PAR)]);
Ex4 = sum((abs(s)).^4);
disp(['E(x4) of result s:' num2str(Ex4)]);
disp(['cve of result s:' num2str(10*log10(cve(s)))]);
disp(['power of result s:' num2str(s'*s)]);
disp('----------------------------------------');


if figswitch
    figure(2);
    start_pow = fft(starts);
    start_pow = conj(start_pow) .* start_pow;
    out_pow = fft(s);
    out_pow = conj(out_pow) .* out_pow;
    subplot(211);hold on;
    if numel(xr_index) > 0
        stem(xr_index, abs(startweights(xr_index)), '--m',...
            'DisplayName', "weights最终值");
    end
    if numel(xi_index) > 0
        stem(xi_index, abs(startweights(xi_index)), 'b',...
            'DisplayName', "信号子载波位置频谱");
    end
    plot(abs(weights), 'r', 'DisplayName', "信号子载波位置频谱");
    legend();
%     legend("weights起始值", "weights最终值", "信号子载波位置频谱")
    title("频谱");
    subplot(212);hold on;
    plot(start_pow, 'r');
    plot(out_pow, '-.b');
    legend("初始点功率谱", "最终功率谱")
    axis([0 128 0 0.1])
    carrierindex = find(1-mask);
    if(isempty(carrierindex))
        titlestr = ['优化前后频谱, carrier:',...
        num2str(0),'-', num2str(0),...
        ', os16Finalcve:', num2str(10*log10(os16Finalcve)),...
        'dB'];
    else
        titlestr = ['优化前后频谱, carrier:',...
        num2str(carrierindex(1)),'-', num2str(carrierindex(end)),...
        ', os16Finalcve:', num2str(10*log10(os16Finalcve)),...
        'dB'];
    end
    title(titlestr);
    
    
    figure(1);
    subplot(211);hold on;
    ylim([0 7e-4]);
    plot(abs(os16StartS), 'b');
    plot(abs(os16OutS), 'r');
    title("最终输出信号")
%     axis([0 2048 0 10e-4]);
    subplot(212);hold on;
    yyaxis left;
    plot(abs(weights(1:N)), 'r');
    plot(abs(startweights(1:N)), 'b');
    title("最终输出的调制系数")
%     xlim(gca, [0 128*3])
    ylim(gca, [0 0.1])
    grid on;
    
    figure(21)
    hold on;
    plot(abs(acf_start), '-r','LineWidth', 1);
    plot(abs((acf_out)), '-.b','LineWidth', 1.5);
%     xlim(gca, [(Os)-200 2048+200])
    legend("起始的ACF", "最终输出信号的ACF");
    title("输出波形的自相关")
    
%     figure(22)
%     histogram(dels_pow(2:end),256);
%     title("子载波功率变化的直方图(仅自由变化部分)")
    
    figure(23)
    hold on
    plot(10*log10(cvesrocord));
    title("cve changes with iterations")

    
%     figure(24)
%     histogram(dels_ampspec(2:end),256);
%     title("频谱幅度起伏的直方图(仅自由变化部分)")
end





% for i=3:Os*N
%     figure(51)
%     
%     delta2 = srecord(i,:) ./ srecord(i-1,:);
%     delta1 = srecord(i-1,:) ./ srecord(i-2,:);
% % %     subplot(211)
% % %     title("s幅值变化");
% % %     plot(abs(srecord(i,:)));   
% % %     ylim([2e-3 6e-3]);
% %     
% %     subplot(211)
% %     title("恒模向量角度随迭代变化趋势");
% %     plot(unwrap(bangles(:,i)));
% %     ylim([-0.15 0.15]);
%         
%     subplot(212)
% %     yyaxis left;hold off;
% %     plot(abs(delta2), 'r-');
% %     hold on;
% %     plot(abs(delta1), 'b-');
% %     plot(abs(delta1.^2), 'b--');
% %     ylim([0.5 1.5]);
%     
%     yyaxis right;hold off;
%     plot(angle(delta2), 'r-');
%     hold on;
% %     plot(2*angle(delta1), 'b-');
%     plot(angle(delta1), 'b--');
%     ylim([-pi/4 pi/4]);
%     yticks([-pi/4 -pi/8 0 pi/8 pi/4]);
%     
%     disp(['第' num2str(i) '次迭代']);
%     
% %     figure(52)
% %     hold off;
% %     plot(abs(srecord(i,:)), 'r');
% %     hold on;
% %     plot(abs(srecord(i-2,:) .* (delta1.^2)), 'g');
% %     plot(abs(srecord(i-2,:)), 'b');
%     
% end



% 研究第一次迭代与最后一次迭代的关系
% scalrefiteb = refiteb * sqrt(N / (N - infnum)) / Os / N;
% scalfirstiteb = firstiteb * sqrt(N / (N - infnum)) / Os / N;
% scalseciteb = seciteb * sqrt(N / (N - infnum)) / Os / N;
% scalstartb = startb * sqrt(N / (N - infnum)) / Os / N;
% 
% deltabbetweenstartandfirst = scalfirstiteb - scalstartb;
% deltabbetweenstartandsec = scalseciteb - scalstartb;
% deltabbetweenstartandlast = scalb - scalstartb;
% 
% figure(47)
% 
% subplot(211)
% hold on;
% 
% deltabbetween1and2 = scalfirstiteb - scalrefiteb;
% deltabbetween2and3 = scalseciteb - scalrefiteb;
% diffofdelta = deltabbetween1and2' * deltabbetween2and3 / norm(deltabbetween1and2, 2) / norm(deltabbetween2and3, 2);
% disp("diff between iteration delta :" + num2str(diffofdelta) + ", norm(diff,2):" + num2str(norm(diffofdelta, 2)));
% 
% 
% 
% plot(abs(fft(deltabbetween1and2)), 'r');
% plot(abs(fft(deltabbetween2and3)), 'g');
% corr = deltabbetween1and2' * deltabbetween2and3 /( norm(deltabbetween1and2, 2) ^ 2);
% plot(abs(fft(15*deltabbetween1and2)), '-.b');
% plot(abs(fft(scalrefiteb)), 'k');
% 
% subplot(212)
% hold on
% plot(abs(scaldeltacoff*fftdeltab(outx_index)), 'r');
% plot(abs(fftb(outx_index)), 'g');

figure(44);hold on;
meanpower = zeros(size(weights));
thpower = zeros(size(weights));

meanpower(:,:) = 0.0945;
thpower(:,:) = 3*0.0945;
plot(abs(weights), 'r', 'LineWidth', 0.75);
plot(meanpower, '--b', 'LineWidth', 1.25);
plot(thpower, '--m', 'LineWidth', 2);
xlabel("index of subcarrier", 'Units', 'centimeters');
ylabel("scale normalized by E(|{\bf\it\fontname{Times New Roman}x}_{\it\fontname{Times New Roman}i}|)");
xlim([0 120])
ylim([0 1.2*thpower(1)])
yticks(0.0945 * [0 0.2 0.4 0.6 0.8 1.0 1.2 1.4]);
yticklabels({'0.0', '0.2', '0.4', '0.6', '0.8', '1.0', '1.2', '1.4'})

set(gca, 'XTickLabel', []);
set(gca, 'YTickLabel', []);
set(gcf, 'unit', 'centimeters', 'position', [10 10 8 6.0]);
set(gca, 'unit', 'normalized', 'position', [0.1200 0.10 0.8650 0.8800]);
set(gca, 'Fontsize',8);
set(gca, 'LineWidth',1,'GridLineStyle','--','GridColor','k','GridAlpha',0.5);

legends = {'modulus of weights', 'E(|{\bf\it\fontname{Times New Roman}x}_{\it\fontname{Times New Roman}i}|)', '\alphaE(|{\bf\it\fontname{Times New Roman}x}_{\it\fontname{Times New Roman}i}|),\alpha=1.2'};
legendhanlder = legend(legends);
set(legendhanlder, 'Box', 'on');
grid on
yticklabels({'0.0', '0.2', '0.4', '0.6', '0.8', '1.0', '1.2', '1.4'})

figure(45);
xlim([-0.1*pi 0.1*pi]);
phi = wrapToPi(angle(h));
histogram(phi, 32);
title("优化前后恒模向量各元素旋转角度分布")


figure(56)
hold on
for i =1:Os*N
    plot(unwrap(deltaangewithchis(i,:)));
end

figure(57)
hold on
histogram(wrapToPi(del_pha), 32);


%% 图片移动

% 对比优化最终结果与其对应恒包络信号频谱间的差异
if figswitch
    figure(2)
    movegui("east");
    
%     figure(45);
%     movegui("north");
    
    figure(1);
    movegui("south");
    
    figure(21);
    movegui("west");
end
figure(54)
% clear all;

%% 幅度归一化研究

% cirmatrix = zeros(Os*N, Os*N);
% for i = 1:Os*N
%     cirmatrix(i, :) = circshift(outs, (1-i));
% end
% d = dftmtx(Os*N);
% cd = conj(d);
% test = cd * cirmatrix * d;
% test = round(test, -8);

%% 加窗方法抑制，经过证明是无效的
% s = outs;
% figure(1);hold on
% plot(abs(s), 'r');
% Ncen = maxpos;
% Nlow = 1683;
% suprat = zeros(Os*N, 1);
% 
% 
% maxr = 0.8;
% lpha = exp (1i * 257 / 2048 * (-10:10)');
% convham = hamming(21) .* lpha;
% 
% hamm = hamming(2*Ncen - 2*Nlow + 1);
% suprat(Nlow: 2*Ncen - Nlow) = - maxr*hamm;
% sups = s + suprat.*flucomplex;
% 
% rweights = fft(sups);
% rweights = rweights(17:128);
% weights = remapi*iweights + remapr*rweights;
% s = ifft(weights, Os*N);
% s = pwr * s / norm(s, 2);
% 
% figure(1)
% plot(abs(s), '--g')
% yyaxis right;
% plot(suprat, '-b')
% yyaxis left;
% plot(abs(sups), '-k')
% 
% figure;hold on;
% plot(abs(fft(sups)), 'r');
% plot(abs(fft(outs)), 'b');
% plot(abs(fft(s - flucomplex)), 'g');
% plot(abs(fft(flucomplex)), '--k');
% 
% figure;
% plot(unwrap(angle(outs)));

%% 设计抑制函数方法

%时域期望的抑制函数
% s = outs;
% supressF = avr ./ abs(s);
% figure(5);hold on
% yyaxis left;
% plot(supressF, '-r');
% fftsuppressF = fft(supressF);
% yyaxis right;
% plot(abs(fftsuppressF), '-b');
% fftsuppressF(65:1984) = 0;
% supressF1 = ifft(fftsuppressF);
% yyaxis left;
% plot(abs(supressF1));



%% 去除特定起伏 
% for i = 1:rounds
%     
%     b = Ad * b + norm(X * b, 2) * c;
%     b = b ./ abs(b);
%     
%     rweights = X * b / norm(X * b, 2);
%     weights = remapi*iweights + remapr*rweights;
%     s = ifft(weights, Os*N);
%     s = sqrt(Os * N) / norm(s, 2) * s;
%     figure(5)
%     plot(abs(s))
% end
%% 非鞍点
% weights = amp .* exp(1i * pha);
% iweights = mapi * weights;
% % load ('iweight_1_16.mat', 'iweights');
% rweights = mapr * weights;
% 
% % 信号分量，全功率，可自由调整信号功率
% powr = sqrt(rweights'*rweights);
% iweights = iweights / powr;
% rweights = rweights / powr;
% weights = remapi*iweights + remapr*rweights;
% s = ifft(weights, Os*N);
% 
% % figure(1);
% % hold on;
% % plot(abs(s));
% 
% X = wr';
% Ad = X' * X;
% % Ad = Ad - diag(diag(Ad));
% 
% % b = wr*rweights;
% % b = b ./ abs(b);
% c = wi * iweights;
% 
% b = wr*rweights + c;
% b = b ./ abs(b);
% 
% % b = c ./ abs(c);
% % bpha = angle(b);
% % bpha = bpha + 0.2*pi*rand(Os*N, 1);
% % b = exp(1j* bpha);
% 
% 
% for i = 1:rounds
%     
% %     rweights = X*b / norm(X*b, 2);
% %     weights = remapi*iweights + remapr*rweights;
% %     s = ifft(weights, Os*N);
% %     delta(i, 1) = cve(s);
%     
%     orgb = b;
%     b = Ad * b + norm(X * b, 2) * c;
%     b = b ./ abs(b);
%     
%     delta(i, 1) = norm(X*b, 2) + real(c'*b);
% %     delta(i, 1) = norm(imag(conj(b) .* (X'*X*b + norm(X*b, 2) * c)) / norm(X*b, 2), 1);
%     
% 
% end
% 
% figure(3);hold on
% % plot(delta(10: rounds), '--b');
% % plot(delta1(40: rounds), 'k');
% % plot(delta2(40: rounds), 'b');
% % plot(delta1(40: rounds) + delta2(40: rounds), 'g');
% 
% rweights = X*b / norm(X*b, 2);
% weights = remapi*iweights + remapr*rweights;
% s = ifft(weights, Os*N);
% 
% figure(2)
% plot(abs(weights), '--b');
% 
% figure(1);
% hold on;
% plot(abs(s), '--k');
% % 
% % figure(2);
% % plot(abs(weights), 'r');
% % 
% disp(['final cve of probable saddle is ', num2str(cve(s))]);
% figure(3)
% % figure(4)

% % 清除外部调用痕迹
clear pca_l1_calloutside

